South Fork Example

Add a watershed outlet point

#Create a point near the stream using lat, long. 
sheds <- tibble(site = c('South Fork'),
                long = c(40.5475185),
                lat = c(-105.6091385)) %>%
  st_as_sf(., coords = c('lat','long'), crs = 4263) %>%
  st_transform(2163)

#Visualize the points
mapview(sheds)

#Little code snippet to make the data folder if you don't have it. 
if(!file.exists('data')){
  dir.create('data')
}

st_write(sheds, 'data/pourpoints.shp', append = F)

## Deleting layer `pourpoints' using driver `ESRI Shapefile'
## Writing layer `pourpoints' to data source 
##   `data/pourpoints.shp' using driver `ESRI Shapefile'
## Writing 1 features with 1 fields and geometry type Point.

Get elevation data from around that area.

# Use elevatr to download data and convert to terra object
elev <- get_elev_raster(sheds, z = 10) %>%
  rast(.)

## Mosaicing & Projecting

## Note: Elevation units are in meters.

# Various ways to plot

mapview(elev) + 
  mapview(sheds)

Save a local DEM Version

#Save local version for whitebox to use
writeRaster(elev, 'data/elev.tif', overwrite = T)

Get hillshade for display purposes

#Get a shillshade map from wbt
wbt_hillshade(dem = "data/elev.tif",
              output = "data/hillshade.tif",
              azimuth = 115)

#Read in the hillshade
hillshade <- rast('data/hillshade.tif')

#Visualize
tm_shape(hillshade)+
  tm_raster(style = "cont",
            palette = "-Greys", 
            legend.show = FALSE)+
  tm_scale_bar()

Hydrologically condition DEM

JP Gannon does a great job on his hydroinformatics website explaining why we need to hydrologically condition a DEM before we can conduct watershed analyses. Basically, DEMs aren’t perfect and they can have artifacts in them that make it so that water doesn’t properly flow downhill. To force water to move the way we think it should through the landscape we can “condition” these DEMs using the functions below.

# Breach depressions (force water to move through ponds/lakes/etc...)
wbt_breach_depressions_least_cost(
  dem = "data/elev.tif",
  output = "data/breached.tif",
  dist = 9,
  fill = TRUE)


# Fill any remaining depressions
wbt_fill_depressions_wang_and_liu(
  dem = "data/breached.tif",
  output = "data/breachfill.tif"
)

Flow Accumulation

Flow accumulation is the key function we use to estimate how much watershed area is draining to a specific point anywhere in the river network. Here we use the D8 flow algorithm, which makes it so 100% of the water is routed to nearby cells. ESRI has a nice article visualizing what this algorithm is doing

#Get flow accumulation
wbt_d8_flow_accumulation(input = "data/breachfill.tif",
                         output = "data/d8fa.tif")
#Get flow direction
wbt_d8_pointer(dem ='data/breachfill.tif',
               output = 'data/d8point.tif')

# read in the data
fa <- rast('data/d8fa.tif') %>%
  log10(.)

#visualize
mapview(fa) + 
  mapview(sheds)

Extract streams with arbitrary 300 cell threshold

In order to make sure that our watershed outlet point, which we arbitrarily and manually extracted lats and longs for, we need to first extract a stream network. We are picking a 300-cell flow accumulation threshold for our ‘stream initiation.’ Because our starting raster is a ~58X58m cell size 300 cells roughly equals 1 km2. This means, we are saying that once a stream has more than 1 km2 of watershed area draining to that point in the stream, we think that it will have a functioning stream channel. You can conduct entire studies to get this threshold right, so this is arbitrary here.

#Extract streams at 1km2 
wbt_extract_streams(flow_accum = "data/d8fa.tif",
                    output = "data/raster_streams.tif",
                    threshold = 300)

#Snap our watershed flowliens to this point. 
wbt_jenson_snap_pour_points(pour_pts = "data/pourpoints.shp",
                            streams = "data/raster_streams.tif",
                            output = "data/snappedpp.shp",
                            snap_dist = 200)

# Read in and check the snap. 
snap_pour <- st_read('data/snappedpp.shp')

## Reading layer `snappedpp' from data source 
##   `/Users/localmatt/Library/CloudStorage/OneDrive-Colostate/ROSS Lab Shared Files/Teaching/523c/Terrain_Analysis/data/snappedpp.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -474945.1 ymin: -478425.3 xmax: -474945.1 ymax: -478425.3
## Projected CRS: NAD27_US_National_Atlas_Equal_Area

mapview(fa) + 
  mapview(snap_pour)

Delineate watershed

Finally after all that! We can delineate our watershed

wbt_watershed(d8_pntr = "data/d8point.tif",
              pour_pts = "data/snappedpp.shp",
              output = "data/south_fork.tif")

south_shed <- rast('data/south_fork.tif')

mapview(south_shed) + 
mapview(sheds)

Extract a terrain metric about this watershed

Generate terrain metric

Delineating a watershed allows us to now get watershed terrain metrics for that watershed. Whitebox tools has dozens of terrain tools you can use See More in the Geomorph Section. Here we will generate and extract one critical variable. Topographic Wetness, which indicates areas in a watershed (or a watershed average) wetness. High values indicate areas that can be marshy/wetlandy and generally wet. Low values (like ridge tops) will be drier parts of the landscape.

#First we need slope
wbt_slope(dem = 'data/elev.tif',
          output = 'data/slope.tif',
          units = 'degrees')

# Now we can generate twi
wbt_wetness_index(sca = 'data/d8fa.tif',
                  slope = 'data/slope.tif',
                  output = 'data/twi.tif')

twi <- rast('data/twi.tif')

mapview(twi)

Crop to watershed and extract average value

# Crop TWI to south_shed area

twi_south <- crop(twi, south_shed) 


twi_south_mask <- mask(twi, south_shed)


mapview(twi_south_mask)

global(twi_south_mask, fun = 'mean', na.rm = T)

##         mean
## twi 2.896648

Terrain Analyses in R

Matthew Ross

2024-04-01

Goals

Terrain Analysis Toolkit

Terrain Analysis Ideas